home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / RLaB / toolbox / mdsmax.r < prev    next >
Text File  |  1994-04-25  |  7KB  |  217 lines
  1. # function [x, fmax, nf] = mdsmax(fun, x, stopit, savit)
  2. #MDSMAX  [x, fmax, nf] = MDSMAX(fun, x0, STOPIT, SAVIT) attempts to
  3. #        maximize the function specified by the string fun, using the
  4. #        starting vector x0.  The method of multi-directional search is used.
  5. #        Output arguments:
  6. #               x    = vector yielding largest function value found,
  7. #               fmax = function value at x,
  8. #               nf   = number of function evaluations.
  9. #        The iteration is terminated when either
  10. #               - the relative size of the simplex is <= STOPIT(1)
  11. #                 (default 1e-3),
  12. #               - STOPIT(2) function evaluations have been performed
  13. #                 (default inf, i.e., no limit), or
  14. #               - a function value equals or exceeds STOPIT(3)
  15. #                 (default inf, i.e., no test on function values).
  16. #        The form of the initial simplex is determined by STOPIT(4):
  17. #           STOPIT(4) = 0: regular simplex (sides of equal length, the default)
  18. #           STOPIT(4) = 1: right-angled simplex.
  19. #        Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
  20. #        If a non-empty fourth parameter string SAVIT is present, then
  21. #        `SAVE SAVIT x fmax nf' is executed after each inner iteration.
  22. #        NB: x0 can be a matrix.  In the output argument, in SAVIT saves,
  23. #            and in function calls, x has the same shape as x0.
  24.  
  25. # This implementation uses 2n elements of storage (two simplices), where x0
  26. # is an n-vector.  It is based on the algorithm statement in [2, sec.3],
  27. # modified so as to halve the storage (with a slight loss in readability).
  28.  
  29. # References:
  30. # [1] V.J. Torczon, Multi-directional search: A direct search algorithm for
  31. #     parallel machines, Ph.D. Thesis, Rice University, Houston, Texas, 1989.
  32. # [2] V.J. Torczon, On the convergence of the multi-directional search
  33. #     algorithm, SIAM J. Optimization, 1 (1991), pp. 123-145.
  34. # [3] N.J. Higham, Optimization by direct search in matrix computations,
  35. #     Numerical Analysis Report No. 197, University of Manchester, UK, 1991;
  36. #     to appear in SIAM J. Matrix Anal. Appl, 14 (2), April 1993.
  37.  
  38. # By Nick Higham, Department of Mathematics, University of Manchester, UK.
  39. #                 na.nhigham@na-net.ornl.gov
  40. # July 27, 1991.
  41.  
  42. #
  43. # Translated to RLaB, Ian Searle
  44. # Febuary 1994.
  45. #
  46.  
  47. # Dependencies
  48. rfile rem
  49.  
  50. mdsmax = function (fun, X, stopit, savit)
  51. {
  52.   local (mu, scale, fmax_old, T, fmax, V, trace, tol, ...
  53.          theta, replaced, size_simplex, alpha, flag_break, ...
  54.          v1, nf, f, ft, msg, x, x0, j, k, m, n, msize)
  55.  
  56.   x = X;    # Copy input
  57.   n = prod(size(x));
  58.   x0 = x[:];    # Work with column vector internally.
  59.  
  60.   mu = 2;      # Expansion factor.
  61.   theta = 0.5; # Contraction factor.
  62.  
  63.   # Set up convergence parameters etc.
  64.   if (!exist(stopit)) { stopit[1] = 1e-3; }
  65.   tol = stopit[1];  # Tolerance for cgce test based on relative size of simplex.
  66.   if (max(size(stopit)) == 1) { stopit[2] = inf(); }    # Max no. of f-evaluations.
  67.   if (max(size(stopit)) == 2) { stopit[3] = inf(); }    # Default target for f-values.
  68.   if (max(size(stopit)) == 3) { stopit[4] = 0; }    # Default initial simplex.
  69.   if (max(size(stopit)) == 4) { stopit[5] = 1; }    # Default: show progress.
  70.   trace  = stopit[5];
  71.   if (!exist(savit)) { savit = []; }            # File name for snapshots.
  72.  
  73.   V = [zeros(n,1), eye(n,n)]; T = V;
  74.   f = zeros(n+1,1); ft = f;
  75.   V[;1] = x0; 
  76.   x = reshape (x0, x.nr, x.nc);
  77.   f[1] = fun (x);
  78.   fmax_old = f[1];
  79.  
  80.   if (trace) { printf("f(x0) = %9.4e\n", f[1]); }
  81.  
  82.   k = 0; m = 0;
  83.  
  84.   # Set up initial simplex.
  85.   scale = max([norm(x0,"i"),1]);
  86.   if (stopit[4] == 0)
  87.   {
  88.     # Regular simplex - all edges have same length.
  89.     # Generated from construction given in reference [18, pp. 80-81] of [1].
  90.     alpha = scale / (n*sqrt(2)) * [ sqrt(n+1)-1+n, sqrt(n+1)-1 ];
  91.     V[;2:n+1] = (x0 + alpha[2]*ones(n,1)) * ones(1,n);
  92.     for (j in 2:n+1)
  93.     {
  94.       V[j-1;j] = x0[j-1] + alpha[1];
  95.       x = reshape (V[;j], x.nr, x.nc);
  96.       f[j] = fun (x);
  97.     }
  98.   else
  99.     # Right-angled simplex based on co-ordinate axes.
  100.     alpha = scale*ones(n+1,1);
  101.     for (j in 2:n+1)
  102.     {
  103.       V[;j] = x0 + alpha[j]*V[;j];
  104.       x = reshape (V[;j], x.nr, x.nc);
  105.       f[j] = fun (x);
  106.     }
  107.   }
  108.   nf = n+1;
  109.   msize = 0;        # Integer that keeps track of expansions/contractions.
  110.   flag_break = 0;    # Flag which becomes true when ready to quit outer loop.
  111.  
  112.   while (1)    ###### Outer loop.
  113.   {
  114.     k = k+1;
  115.  
  116.     # Find a new best vertex  x  and function value  fmax = f(x).
  117.     fmax = max (f); j = maxi (f);
  118.     V[;1,j] = V[;j,1]; 
  119.     v1 = V[;1];
  120.     if (!isempty(savit)) { x = reshape(v1, x,nr, x.nc); write ("savit", x,fmax,nf); }
  121.     f[1,j] = f[j,1];
  122.     if (trace)
  123.     {
  124.       printf("Iter. %2.0f,  inner = %2.0f,  size = %2.0f,  ", k, m, msize);
  125.       printf("nf = %3.0f,  f = %9.4e  (%2.1f)\n", nf, fmax, ...
  126.              100*(fmax-fmax_old)/(abs(fmax_old)+eps));
  127.     }
  128.     fmax_old = fmax;
  129.  
  130.     # Stopping Test 1 - f reached target value?
  131.     if (fmax >= stopit[3])
  132.     {
  133.       msg = "Exceeded target...quitting\n";
  134.       break  # Quit.
  135.     }
  136.  
  137.     m = 0;
  138.     while (1)   ### Inner repeat loop.
  139.     {
  140.       m = m+1;
  141.   
  142.       # Stopping Test 2 - too many f-evals?
  143.       if (nf >= stopit[2])
  144.       {
  145.         msg = "Max no. of function evaluations exceeded...quitting\n";
  146.         flag_break = 1; 
  147.         break  # Quit.
  148.       }
  149.   
  150.       # Stopping Test 3 - converged?   This is test (4.3) in [1].
  151.       size_simplex = norm(V[;2:n+1] - v1[;ones(1,n)],"1") / max([1, norm(v1,"1")]);
  152.       if (size_simplex <= tol)
  153.       {
  154.         sprintf(msg, "Simplex size %9.4e <= %9.4e...quitting\n", ...
  155.                 size_simplex, tol);
  156.         flag_break = 1; 
  157.         break  # Quit.
  158.       }
  159.   
  160.       for (j in 2:n+1)      # ---Rotation (reflection) step.
  161.       {
  162.         T[;j] = 2*v1 - V[;j];
  163.         x = reshape (T[;j], x.nr, x.nc);
  164.         ft[j] = fun (x);
  165.       }
  166.       nf = nf + n;
  167.   
  168.       replaced = ( max(ft[2:n+1]) > fmax );
  169.   
  170.       if (replaced)
  171.       {
  172.         for (j in 2:n+1)    # ---Expansion step.
  173.         {
  174.           V[;j] = (1-mu)*v1 + mu*T[;j];
  175.           x = reshape (V[;j], x.nr, x.nc);
  176.           f[j] = fun (x);
  177.         }
  178.         nf = nf + n;
  179.         # Accept expansion or rotation?
  180.         if (max(ft[2:n+1]) > max(f[2:n+1]))
  181.         {
  182.           V[;2:n+1] = T[;2:n+1];
  183.           f[2:n+1] = ft[2:n+1];    # Accept rotation.
  184.         else
  185.           msize = msize + 1;    # Accept expansion (f and V already set).
  186.         }
  187.       else
  188.         for (j in 2:n+1)    # ---Contraction step.
  189.         {
  190.           V[;j] = (1+theta)*v1 - theta*T[;j];
  191.           x = reshape (V[;j], x.nr, x.nc);
  192.           f[j] = fun (x);
  193.         }
  194.         nf = nf + n;
  195.         replaced = ( max(f[2:n+1]) > fmax );
  196.         # Accept contraction (f and V already set).
  197.         msize = msize - 1;
  198.       }
  199.   
  200.       if (replaced) { break }
  201.       if (trace && rem(m,10) == 0) 
  202.       { 
  203.         printf("        ...inner = %2.0f...\n",m);
  204.       }
  205.     }        ### Of inner repeat loop.
  206.  
  207.     if (flag_break) { break }
  208.  
  209.   }    ###### Of outer loop.
  210.  
  211.   # Finished.
  212.   if (trace) { printf(msg); }
  213.   x = reshape (v1, x.nr, x.nc);
  214.  
  215.   return << x = x; fmax = fmax; nf = nf>>;
  216. };
  217.